Halide Adsorption on Single-crystal Silver Substrates: 

Dynamic Simulations 
and ab-initio Density-functional Theory 

Q ; S. J. Mitchell, Sanwu Wang, and P. A. Rikvold 

Center for Materials Research and Technology, 
School of Computational Science and Information Technology, 

and Department of Physics, 
Florida State University, Tallahassee, Florida, 32316-4120, USA 

February 1, 2008 



o 



o 
o 



. 

£h ■ We investigate the static and dynamic behaviors of a Br adlayer electro chemically deposited 
onto single-crystal Ag(100) using an off-lattice model of the adlayer. Unlike previous studies 
using a lattice-gas model, the off-lattice model allows adparticles to be located at any position 
t— ( \ within a two-dimensional approximation to the substrate. Interactions with the substrate are 
^ [ approximated by a corrugation potential. Using Density Functional Theory (DFT) to calculate 
■ surface binding energies, a sinusoidal approximation to the corrugation potential is constructed. 
A variety of techniques, including Monte Carlo and Langevin simulations, are used to study the 
behavior of the adlayer. The lateral root-mean-square (rms) deviation of the adparticles from 
the binding sites is presented along with equilibrium coverage isotherms, and the thermally 
activated Arrhenius barrier-hopping model used in previous dynamic Monte Carlo simulations 
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1 Introduction 

o : 

\ Electrochemical adsorption processes in monolayer or submonolayer systems contain essentially 
| four types of microscopic processes: lateral surface diffusion, adsorption/desorption, molecular 
^ ■ reorientation, and reaction or charge transfer. Often, the free-energy landscape for these sys- 
■ terns has a very large number of local minima, each of which corresponds to a well-localized 
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state in the full phase space. For example, adparticles are often well localized on the crystal 
substrate, spending the majority of their time near specific adsorption sites. In these systems, 
the time scale for transitions between neighboring free-energy minima is associated with hop- 
ping across a free-energy barrier that corresponds to a saddle point in the free-energy landscape. 
Dynamic Monte Carlo simulations exploit these well-localized states and simulate the time evo- 
lution of the system as a series of thermally activated stochastic transitions over the free-energy 
saddle points. 

In this paper we investigate the validity of the localized-state approximation and the com- 
monly used Arrhenius form for the stochastic transition rate between local free-energy minima 
by constructing an off-lattice model for the adsorption of Br onto the (100) surface of an Ag 
single crystal. This system has the advantage of simplicity, in that the only microscopic pro- 
cesses are lateral surface diffusion and adsorption/desorption. In this system, the only reaction 
or charge-transfer process is associated with adsorption/desorption, and reorientation is not 
important for atomic adsorbates. 

The equilibrium properties of Br adsorbed onto Ag(100) have been extensively studied both 
by experiments [0, |^, |3[] and by theory ||, [|, || [| . The time evolution of a lattice-gas approxi- 
mation to the adsorbate dynamics has also been extensively studied by dynamic Monte Carlo 
simulations [|, 0, |J- The lattice-gas approximation used in Refs. ||, ||, |], [7|, §] explicitly assumes 
well-localized surface states, which is not an unreasonable assumption since this system is known 
to display only commensurate phases 0. The present paper is primarily concerned with test- 
ing the well-localized state and Arrhenius transition-rate assumptions of Dynamic Monte Carlo 
for surface diffusion processes. However, the computational techniques used here can also be 
applied to adsorption/desorption and molecular reorientation barrier- hopping processes. The 
off-lattice model developed here, which does not assume well- localized surface adsorption sites, 
can also be used for systems which do not have strongly localized adsorption sites, such as 
the incommensurate phases observed for Br/Au(100) f9j. Some preliminary results of the work 
presented here were included in Ref. [[TH|]. 

The remainder of this paper is organized as follows. In Sec. |2|, Density Functional Theory 
(DFT) pseudopotential calculations are used to calculate the surface binding energies and 
deduce an approximate corrugation potential for Br on Ag(100). In Sec. |3], an off-lattice model 
is constructed using the corrugation potential, and the well-localized surface-state assumption 
is tested by equilibrium simulation of this off-lattice model. In Sec. |J the thermally activated 
stochastic Arrhenius transition-rate assumption is tested for diffusion processes by simulating 
single-particle trajectories within the corrugation potential by Langevin simulations. In Sec. [5|, 
our conclusions are presented, and in the Appendix, the grand-canonical off-lattice Monte Carlo 
method is explained in detail. 
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2 Corrugation Potential 



In situ Surface X-ray Scattering (SXS) [§] and X-ray Absorption Fine Structure (XAFS) || 
measurements of Br deposited on Ag(100) in solution, as well as Low Energy Electron Diffrac- 
tion (LEED) measurements Jll]| in vacuum confirm that Br adsorbs at the four-fold hollow sites 
of the substrate, indicating that the corrugation potential has its minima at these points, both 
in vacuum and in solution. However, knowledge of the binding site alone is not enough to de- 
termine the symmetry or amplitude of the approximate corrugation potential, and therefore ab 
initio Density Functional Theory (DFT) pseudopotential calculations were used to determine 
the parameters of the corrugation potential. 

We point out that one must be careful in selecting surface models for ab initio calculations, 
as the binding energy is very sensitive to the way the surface is modeled (cluster or slab models). 
This point is explained in further detail in Ref. We emphasize that calculations with small 
metal clusters should be used with caution since the electronic structures of a small metal 
cluster and a surface of the same metal are usually very different. Small cluster models are 
useful in providing detailed information about the local bonding. However, convergence of the 



binding energy with the cluster size is usually difficult to achieve fll3"|, and cluster calculations 



can even lead to predictions in direct conflict with experimental results, such as for the binding 
site of Br on Au(100) []14[ . In order to simulate metal surfaces accurately large clusters or 
extended surface models (e.g., slab models) are therefore needed. In contrast to metals, small 
clusters are usually quite good models for semiconductors since the electronic states and the 



chemical bonds are much more localized in this case |L5 . 

In our DFT calculations, the Ag(100) surface was modeled by repeated slabs with seven 
metal layers of nine Ag atoms separated by a vacuum region equivalent to seven metal layers, 
corresponding to a 3 x 3 surface cell. Br was placed symmetrically on both sides of the slab. 
Total energies for the hollow, bridge, and on-top configurations, as well as for configurations 
with Br at the midpoints between the hollow and bridge, the hollow and on-top, and the bridge 
and on-top sites were calculated. All the metal atoms were located at their bulk positions, with 
the equilibrium bulk lattice constant of 4.17 A determined by our calculations, as compared 



to the experimental value of 4.09 A |16| . Test calculations showed that relaxation of the Ag 
positions had a negligible effect on the total-energy differences between the hollow and bridge 
configurations, and we therefore relaxed only the Br position in the direction perpendicular to 
the surface. 

Our ab initio total-energy calculations were performed within density functional theory with 
the Vienna ab-initio simulation package (VASP) |T7], |Tg| , using the pseudopotential method and 
a plane- wave basis set. The exchange-correlation effects were treated with generalized gradient- 
corrected exchange-correlation functionals (GGA) in the form of Perdew and Wang [TT3, |21J, and 
we adopted the Vanderbilt ultrasoft pseudopotentials generated by Kresse and Hafner |^T], [22| . 
The calculations were conducted with a plane-wave energy cutoff of 20 Ry and 9 special k 
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points in the irreducible part of the two-dimensional Brillouin zone. The convergence of the 
total-energy differences (i.e., the binding-energy differences) between different configurations 
was checked with 9 and 16 special k points using supercells containing five metal layers. We 
found that the total-energy differences obtained from calculations with 9 and 16 special k points 
were within a few meV. Extended convergence checks for Br/Ag(100)-c(2 x 2) surfaces suggested 
that the use of a cut-off energy of 20 Ry and supercells with a seven-layer Ag slab separated by 
a vacuum gap equivalent to seven Ag layers resulted in binding-energy differences with errors 
of only a few meV |12| . We therefore believe that the total-energy differences reported in this 
paper were within 10 meV of the values that would be obtained with higher cutoff energies, 
more k points, or thicker supercells. 

The results obtained from the calculations are summarized in Table |T|. The four-fold hollow 
configuration is found to be the most stable one. This is in agreement with the experimental 
observations. The total energy of the bridge configuration is higher by 154 meV than the hollow 
configuration. The on-top configuration is the least stable, with a total energy significantly 
higher than both the hollow and bridge configurations. 

Motivated by the results of these DFT supercell calculations, we approximated the Ag(100) 
surface by a two-dimensional sinusoidal corrugation potential containing L x L four-fold hollow 
adsorption sites, 



U(x,y) = — 



2ttx\ (2ny\ 



cos [ | + cos 
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where x G (0, aL] and y G (0, aL] are the positions on the surface in the two perpendicular 
lattice directions, a = 2.889 A is the two-dimensional lattice constant for Ag(100), and 
A = -Ebridge — -^hollow = -^top — -^bridge gives the amplitude of the corrugation potential. This 
sinusoidal corrugation potential contains the expected rotational symmetries, but it also has 
up/down symmetry, which the actual corrugation potential does not have. We expect this 
additional symmetry to have very little effect on the simulations. 

The differences between the actual DFT-calculated energies and the approximation, Eq. ([[]), 
with A = 150 meV are shown in the last column of Table |I[ The value of A w 150 meV is 
suggested by bridge — -^hollow, as shown in the third column of Table |l|. The difference between 
the two potentials is more pronounced than expected from our earlier c(2 x 2) calculations 



121 . The potential well around the four-fold hollow site as obtained by the DFT calculations is 
much narrower than that given by the sinusoidal form, and this difference could be important. 
However, in this work we are primarily interested in demonstrating the general simulation meth- 
ods required for an off-lattice simulation, for which the sinusoidal approximation is perfectly 
adequate. 

Figure |TJ shows a grayscale plot of U(x,y). The dynamic lattice-gas model presented in 
Refs. [7], [8| implicitly assumed a corrugation potential of this form, with A = 100 meV as 
the nearest-neighbor diffusion barrier. This somewhat lower amplitude was motivated by the 



earlier ab initio cluster calculations of Ref. JT^]. As shown in Table [T], our more recent DFT 
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supercell calculations suggest A « 150 meV 
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3 Equilibrium Off-lattice Simulations 

In this section we test the lattice-gas approximation of the Br adlayer previously assumed in 
Refs. [|5|, H |7|, M. Since lattice-gas models assume well-localized surface states and therefore have 
only discrete adparticle positions, they are incapable of reproducing behaviors which depend on 
lateral degrees of freedom of the adparticle within the adsorption well. We therefore present an 
off-lattice model for the Br adlayer using the approximate sinusoidal corrugation potential based 
on the adsorption energies calculated by DFT in Sec. |2|. In this off-lattice model, adparticle 
positions have continuous degrees of freedom on a two-dimensional surface. 

3.1 Model 

The off-lattice model consists of two parallel planes separated by an unspecified distance, as 
shown in Fig. |^. One plane represents the effective Ag(100) surface (the surface layer), which 
includes a corrugation potential for the binding of Br to the surface. The other plane represents 
the entire volume of solution near the Ag(100) surface, from which a Br ion can adsorb within 
one algorithmic step (the solution layer). The number of Br ions in the solution layer is iV BO i, and 
the number of adsorbed Br is iV sm .f. Both square planes have area (La) 2 . As with the lattice-gas 
models, periodic boundary conditions are used to eliminate edge effects in the relatively small 
systems it is possible to simulate. 

Each Br particle has a position (x, y) within one of the layers. Diffusion moves alter the (x, y) 
positions of the particles within a layer, while only adsorption/desorption moves can transfer 
particles between the layers. The grand-canonical Hamiltonian for this model is 

1 -V sur f 

W = -- fa + J2 U Vi) ~ Wurf , (2) 
1 i+j i=l 

where i and j index adsorbed Br, | YUj=j i s a sum over an Br pairs where the factor of 1/2 is 
included to prevent double-counting, faj is the lateral Br-Br interaction energy for the pair i,j 
within the adlayer, J2i=L f is a sum over all adsorbed Br particles, U(xi,yi) is the value of the 
corrugation potential for the z-th Br, and /2 is the electrochemical potential. All energies are 
measured in units of meV/particle or meV/pair, and all positions and distances are measured 
in units of A. To simplify the notation, all energy units are simply denoted as meV. The sign 
convention used here is such that negative implies repulsion, more negative U represents 
stronger binding energy, and p, > favors adsorption. 

Motivated by the success of the lattice-gas models || 0) 0? we represent the lateral 
interactions for large separations as a 1/r 3 repulsion and for separations less than the ionic 
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diameter as a truncated Lennard- Jones potential. The 1/r 3 potential is most likely associated 
with both dipole-dipole and lattice-mediated interactions. Br in the solution layer are non- 
interacting. The lateral interaction potential is 

-0nnn(ir) 3 T < Aon 

(f) 3 " Aon<r<5a > ( 3 ) 

r > 5a 




'nnn ^ 5 J 



where r is the separation between the members of an interacting adsorbed Br pair, Aon is 
the ionic diameter of Br (Aon — 2i? ion where i? ion = 1.94 A fT5|]), </> nnn = —26 meV (as 
determined by fitting lattice-gas simulations to experimental adsorption isotherms in Ref. 0) 
is the value of <p a ^ next-nearest neighbor distances (r = ay/2). The short-range Lennard- 
Jones parameters, a = i? ion and e = —(p nnn (a\/2/R ion ) 3 /S, are found by requiring that cf>(r) be 
continuous everywhere and that d0(r)/dr be continuous at r = Aon- Equation @ is plotted 
in Fig. |3|. 

In the weak-solution approximation, the electrochemical potential p, is related to the tem- 
perature, the concentration of Br ions in solution, and the electrode potential by 

C 

Ji = fio + k B T In — - 7eA (4) 

where JIq is a reference potential, k B is Boltzmann's constant, T is the absolute temperature 
(ksT = 25 meV at room temperature), C is the concentration of Br ions in solution, Co 
is a reference concentration of Br ions in solution, 7 is the electrosorption valency, e is the 
elementary charge unit, and E is the electrode potential in mV. Since the off-lattice models 
considered here involve a form of the solution concentration, iV BO i, all isotherms require a shift 
of Jl for comparison with the lattice-gas isotherms, as further explained in the Appendix. 

3.2 Methods 

The equilibrium properties of a system are independent of the dynamical processes. Thus, we 
here use the standard Metropolis Monte Carlo (MC) simulation |23| to describe the equilibrium 



properties of the adlayer at room temperature. The essence of the MC method is that, given 
the current state I of the system, a new state F is randomly chosen by making a small change 
to the state I. Once the energies Ej and E F of states / and F are calculated using Eq. (g), the 
new state F is accepted with probability 

P(F\I) = min [l,exp (-(E F - Ej)/k B T)} . (5) 

We discuss two different types of off-lattice MC simulations. The first simulations are per- 
formed within a fixed- coverage (canonical) ensemble (i.e., iV sur f constant). From these simu- 
lations we measure the lateral root-mean-square (rms) displacement of the adparticles from 
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the four-fold hollow sites as a function of the coverage and the corrugation amplitude A. The 
second set of simulations is performed within the grand-canonical ensemble (i.e., fixed p), and 
the coverage isotherms are calculated for different values of the corrugation amplitude A. Since 
the fixed-coverage simulations are simpler, we discuss these here. The grand-canonical MC 



method [Q is discussed in detail in the Appendix. 

The fixed- coverage simulations begin by placing iV sur f = QL 2 adparticles randomly on the 
surface. To simplify the initial MC relaxations, all particles are placed on the same c(2 x 2) 
sublattice. New trial configurations are generated by selecting one of the iV sur f particles at 
random. Call this particle i. A random displacement vector, dr (which is uniformly distributed 
within a circle of radius i?i on /4), is added to r*j, the energy change is calculated, and the trial 
configuration is accepted with probability given by Eq. (|). 

The grand-canonical MC simulations are significantly more complex than the fixed- coverage 
simulations due to the added difficulties of satisfying detailed balance under the adsorption/desorption 
and lateral diffusion processes. These difficulties can be overcome by using a "ghost-particle" 
simulation method like that of Ref . |24j , which includes both the surface plane and the solution 



plane shown in Fig. 0. A full description of the grand-canonical off-lattice MC method is given 
in the Appendix. 

3.3 Results 

Using canonical MC simulations at room temperature, the lateral rms displacement of the Br 
adparticles from the four-fold hollow sites were measured as a function of coverage and A, as 
shown in Fig. ||. Two behaviors are obvious. First, the rms displacement is larger for smaller 
corrugation amplitudes. Second, the rms displacement is smallest for coverages above the crit- 
ical coverage for the c(2 x 2) structure (6 C ~ 0.37 [0]). This second observation is consistent 
with results of in situ X-ray scattering experiments, which indicate that the lateral rms dis- 
placement is approximately 0.25 A at low coverages [25]. Comparison of these experimental 
results with Fig. £| suggest that A « 100 meV. 

Using the grand-canonical simulation method, off-lattice equilibrium-coverage isotherms for 
0nnn = — 26 meV are shown in Fig. |5] for various corrugation amplitudes A, along with the 
corresponding isotherm from the lattice-gas model. In this figure, the low ft resolution makes 
identifying the phase transition between the disordered low-coverage phase and the ordered 
c(2 x 2) phases || [7| rather difficult. The differences are most notable in the disordered phase 
around = 1/4. For large amplitudes of the corrugation potential (A > 150 meV), the 
disordered region is very linear as seen in both experimental isotherms and the lattice-gas 
isotherm J7|. For lower values of A, the isotherm in the disordered region has a pronounced 
curvature. No attempt to vary has been made because of the large computational demands 
such a study would make. Each of the off-lattice isotherms required approximately 6 weeks of 
CPU time on a 0.5 GHz PC, as compared to just a few days of CPU time for the lattice-gas 
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isotherm shown in Fig. |5|. As discussed in Sec. §, DFT calculations suggest A « 150 meV Jlj ; 
however, comparison of the rms lateral displacement with experiment suggests A m 100 meV. 
The apparent discrepancy between these two estimates may be resolved by a more accurate 
approximation for the shape of the corrugation potential. 

Both the lateral rms displacement and the coverage isotherms indicate that for the reasonable 
value of A 150 meV, the surface states are well localized near the four-fold hollow sites. Thus 
the first approximation made in dynamic Monte Carlo, well-localized states, is valid for surface 
diffusion. 



4 Arrhenius Behavior 

In addition to assuming a lattice-gas treatment of the adlayer, the previous dynamic Monte 
Carlo studies of Br/Ag(100) have assumed that transitions between lattice-gas states are well 
described by thermally activated stochastic barrier hopping || |7|, |§ . In particular, the Arrhe- 
nius form for the transition rates was assumed: 

R(F\I) = u eS exp (—A cfi /k B T) , (6) 

where R(F\I) is the probability per unit time of undergoing a transition between a lattice-gas 
state I and a lattice-gas state F, u e s is an effective attempt frequency, and A eff is the effective 
free-energy barrier. In general, A e g has a dependence on the energies of the states I and 
F, a dependence on the bare barrier of the process, which generally is different for diffusion 
and adsorption/desorption, and it also includes a dependence on entropic effects, such as the 
accessibility of the intermediate transition state. 

To investigate the validity of the Arrhenius assumption, we must perform dynamic simula- 
tions which consider both the continuous particle positions and the kinetic energy. Conventional 
dynamic Monte Carlo ignores the kinetic energy. The ideal approach would be a full molecular- 
dynamics treatment of the substrate atoms, adsorbed Br ions, and all of the atoms in the 
solution. However, such a simulation would require not only knowledge of the very complex in- 
teraction potentials, but also an enormous computational effort. Since the system is dominated 
by thermally activated barrier hopping, it is unlikely that conventional molecular-dynamics 
simulations could reach the necessary long time scales. The full simulation might be feasible 
within the Hyperdynamics Method of Voter pH] , although Hyperdynamics assumes Arrhenius 
type behavior. Thus, we here present a simplified model of the dynamics which still includes 
most of the essential behavior, such as kinetic energy (ignored in MC dynamics) and thermal 
noise. However, as with the lattice-gas models, there is no explicit treatment of the solution 
and substrate atoms. 

We investigate the validity of the Arrhenius assumption by performing Langevin simulations 
of a single Br particle on the surface. The behavior of only a single particle is used to avoid 



S 



lateral Br-Br interactions, thus making the calculations completely general for any system with 
a corrugation potential well approximated by Eq. ([]]). We then analyze the single-particle 
trajectory in different parameter regimes to quantify the Arrhenius behavior and deviations 
from it. 



4.1 Langevin Simulations 

To simplify the simulation, we assume a static substrate, where the average positions of the 
Ag atoms give rise to the corrugation potential of Eq. ([!]), and the fluctuations of the Ag 
atoms about their average positions are included by introducing drag and random-noise terms 
which represent collisions with phonons in the Ag crystal. The presence of water is replaced by 
additional drag and random-noise terms, which are caused by collisions with water molecules. 

We only investigate the lateral diffusion of a single adparticle on the surface, and we therefore 
ignore adsorption/desorption processes. The drag and noise terms from both substrate and 



solution are combined into a Langevin equation of motion |27j governing the lateral motion of 
the Br particle within the surface layer 

mif(t) = -yU(x{t)) - av{t) + ff(t) , (7) 

where m is the effective mass of a Br ion, t is time, v(t) is the time derivative of the velocity 
v(t) of the Br adparticle at time t, x(t) is the two-dimensional position of the adparticle at 
time t, —\/U(x(t)) is the force on the adparticle due to the corrugation potential, a is the 
drag coefficient, and ffi(t) is a Gaussian white random force. For simplicity, we assume that m 
is equal to the mass of a Br atom, and for convenience the values of a are reported in units 
of 80 x 10~ 15 kg/s, corresponding to the system of dimensionless units used in the numerical 
simulations. 



The Gaussian random force has the following properties [27] : 



where the subscripts denote the x and y components of the random force, t and t' are two 
arbitrary times, (■) represents an ensemble or time average, <7q is the variance of a Gaussian 
distribution, and S(-) is the Dirac delta function. The width of the Gaussian noise distribution 
is related to the friction coefficient a and the temperature T by the fluctuation-dissipation 



theorem [27], which requires that 



a G = J2k B Ta . (9) 
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Equation (|7|) is solved iteratively by a first-order integration method, 



v(t + At) = v(t) - yU{x(t))At/m - av(t)At/m + ff(t)vAt/m , 



(10) 



where At is the magnitude of the discrete time step, here taken as At = 1.66 x 10~ 15 s. 
The simulations were not numerically stable for At significantly larger than this value. For 
smaller At, no significant differences were seen. The y/At factor in the term ffi(t)v~At/m is a 
commonly known result of stochastic noise integration |27j and essentially corresponds to the 
noise performing a random walk as a function of time. The new positions are given by 



x (t + At) = x(t) + v{t)At 



4.2 Trajectory Results 

Langevin simulations of a single-particle trajectory were performed at room temperature. We 
define the crossing time, r, as the time spent within the current unit cell of the Ag(100) lattice 
before moving into another adsorption well. The particle was placed on the surface, and after 
throwing away the first few crossings, a time sequence of r was measured including 10,000 
crossings. A typical trajectory is shown in Fig. |6], superimposed on the corrugation potential 
and a grid indicating the unit cells. The bridge site is seen to be the dominant location for 
hopping. Though it is difficult to see in the static image of Fig. three behaviors are easily 
seen in dynamic images of the trajectories: 1) Arrhenius-like single, independent hops between 
cells, 2) crossing and immediate recrossing back into the previous cell, and 3) crossing multiple 
bridge sites in a very short time, which we call "running." The dynamic MC methods used in 
the lattice-gas simulations included only Arrhenius-like crossings. 

The Arrhenius behavior and the small-r deviations from it can be seen in the probability 
density of r, shown in Fig. [7[ The Arrhenius behavior is observed as an exponential tail in 
the distribution, which occurs for sufficiently large r, regardless of the value of the damping 
a (a > 0) or the corrugation amplitude A (A > HbT). However, a large number of crossings 
occur in the short-r regime. Furthermore, in the weakly damped limit (very small a), the 
short-r crossings should be dominated by running, while in the strongly damped limit (very 
large a), crossing/recrossing behavior should dominate the small-r crossings. 

Before examining the large-r Arrhenius behavior in greater detail, we quantify the amount 
of the three crossing behaviors by examining the probability to cross in a particular direction, 
relative to the direction of the previous crossing. When the damping is small (the thermalization 
time is long), then a particle which has sufficient kinetic energy to hop one barrier will also have 
sufficient kinetic energy to immediately hop a second barrier. This is the cause of running. In 
contrast, if the thermalization time is short (large damping), then the small-r crossings should 
correspond to recrossing back into the previous cell. Figure. |8| shows the directional crossing 
probabilities for A = 150 meV and various values of a as functions of r. For a = 1, the small-r 
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regime is dominated by running, but for a = 10 and a = 100, the small-r regime is dominated 
by crossing/recrossing events. 

To further elucidate the long-time Arrhenius behavior, we next examine the temperature 
dependence of the crossing rate. Because the exponential behavior occurs only for large cross- 
ing times, we only examine crossing times larger than a minimum cutoff, r c = 8.3 ps. We 
assume the standard Arrhenius form for the crossing rate given in Eq. (Q). Arrhenius be- 
havior assumes a Poisson process with a constant rate which has the probability distribution 
P(t) = Rexp [— R(t — r c )], for r > r c . After integrating tP(t) for r > r c , we obtain the 
expression 

(r) T >r c = r c + R- 1 = t c + v'g 1 exp [A eS /{k B T)} , (12) 

where (r) T > Tc represents the average crossing time such that r > r c . From this expression, both 
the effective barrier and the attempt frequency can be found from a plot of In [(r) r > Tc — r c )] vs 
(fceT) -1 , as shown in Fig |9|. In this plot, the intercept is — ln(z/), and the slope is -Ebarr- A 
summary of the effective Arrhenius parameters is shown in Table || and Table |3|. Significant 
deviations from the assumption made in Refs. |7|, [], § that -Ebarr = A are seen for the weaker 
damping. Without knowledge of the actual value of a, no prediction of v c q can be made; 
however, we note that for all values of a and A studied here, z/ er r is on the order of 10 12 s _1 . 

As a final comment, the repulsive lateral interactions between Br in the adlayer should alter 
the crossing behavior in the presence of other Br particles. However, these effects have yet to 
be studied. 

5 Conclusions 

To investigate the validity of the lattice-gas approximation or well-localized surface-state as- 
sumption, we have constructed an off-lattice model using a corrugation potential approximated 
on the basis of Density Functional Theory (DFT) surface-binding-energy calculations. Equilib- 
rium Monte Carlo (MC) simulations in both the canonical (fixed coverage) and grand-canonical 
(fixed fx) ensembles were performed. Both the rms lateral displacement from the binding sites 
and the coverage isotherms suggest that for the most reasonable value of the corrugation poten- 
tial, A m 150 meV, the surface states are well localized, thus justifying the previous lattice-gas 
models. 

The stochastic Arrhenius transition rate assumption was also tested using Langevin simu- 
lations of single-particle trajectories. These trajectories are general for any system where the 
corrugation potential is well approximated by the sinusoidal form used in this paper. For long 
time-scale transitions, the transition rate was found to be well described by the Arrhenius rate. 
However, crossing/recrossing and running behaviors were observed for short-time transitions, 
suggesting that future dynamic Monte Carlo simulations should include such behaviors. Un- 
fortunately, the value of the phenomenological damping parameter is unknown for Br/Ag(100), 
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and it is therefore unknown whether the short-time transitions for the physical system fall into 
the crossing/recrossing or running regimes. 

Future theoretical, computational, and experimental efforts should concentrate on determin- 
ing the damping parameter as well as further justification of the dynamic Monte Carlo models. 
Finally, we point out that other processes, such as adsorption/desorption and molecular reori- 
entation, can be studied with similar off-lattice models, where the two-dimensional corrugation 
potential is generalized to three dimensions. 
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Appendix 

In this Appendix we discuss the Monte Carlo methods used to simulate the equilibrium proper- 



ties of the off-lattice model. This method is identical to the ghost-particle method of Ref. |24 . 
However, that reference contains several errors, and the correct details of the method are given 
here. 

The grand-canonical partition function is 

oo 

Z= exp (P»N)Z N , (13) 

N=0 

where N is the number of particles in the adlayer, Zq — 1, and 

Z N = h~ Nd (N\)~ l J d d p x ■ ■ ■ d d p N d d n ■ ■ ■ d d f N exp {-(3Hn) , (14) 

where Pi and r*j are the two-dimensional momentum and position, respectively, of the ith par- 
ticle, H.n is the Hamiltonian of the iV-particle system, including both kinetic and potential 
energies, d is the spatial dimension, (3 = 1/ksT, and h is Planck's constant. For this work, 
d = 2, and we make no distinction between /x and /2. 

As iV —>■ oo, the typical distance between neighboring particles goes to zero. In this limit, 
the interaction energy will go to infinity (see Chapter 4), and thus the exponential factor in 
Eq. (|i~4|) will go to zero. Thus, we can consider some maximum number of particles, M, such 
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that the probability of having a configuration with N > M is negligible. This leads to the 
approximation of Eq. (|T3"D, 

M 

"£exp(PiMN)Z N . (15) 

N=0 

However, this is the partition function for a canonical ensemble describing the system consisting 
of the surface layer and the solution together. The remainder of this Appendix will be concerned 
with constructing such a canonical ensemble which correctly approximates a grand-canonical 
ensemble for the surface layer. 

Consider a system of M particles with continuous momentum and position degrees of free- 
dom. We add the additional discrete degree of freedom, c, which is 1 for particles in the surface 
layer and for particles in the solution layer. Particles in the solution layer are sometimes 
referred to as "ghost" particles. This system has 2dM continuous degrees of freedom and M 
discrete degrees of freedom. The number of particles in the adlayer is then 

M 

iv = E c ^ ( 16 ) 
i=i 

where q denotes whether the zth particle is in the solution (0) or the surface (1) layer. The 
Hamiltonian for this system is 

M 2 M-l M M 

Ktot = E f- + E E U ™t{\n - r,\)c tCj + Oe*(ri)cj + $(iV, M, //, /3, V) , (17) 
i=i ^ m i=i j=i+i i=i 

where m is the effective mass of a particle, Ui nt (\fi — fj\) is the interaction potential between 
two particles, t/ ext (r;) is an external potential, here equivalent to the corrugation potential, and 
$(iV, M, pi, (3, V) is a correction factor which we will now derive. Note also that particles in the 
solution layer have only kinetic energy. 
The partition function for this system is 

Ztot = h~ Md (M\y 1 E E ■ • • E / dVi • • • d d p M d d n ■ ■ ■ d d f M exp (-(3H tot ) . (18) 

ci=0c2=0 cm=0 

The correction factor in Eq. (|17]) can be found by equating Eq. ( |18D and Eq. (pT5|) . We first 
rewrite Eq. 



M 
N 

(19) 



Ztot= h- Md {M\)- l Y. h LA AT )Jd d p 1 ---d d p M d d f 1 ---d d r M 



xexp (-13 [Efit H + Eil! 1 £f =i+1 U int (\n - r^c 
+ Eli U ext (n)ci + ${N, M, fi, (3, V) ' 



which can be rewritten in this way for the following reasons: 
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1. The particles can always be renumbered such that the first N are surface-layer particles. 
This means it is superfluous to write Q,Cj in the surface sums, since they are unity by 
definition. 



2. The number of ways to assign the occupation variables q for N surface particles is 



M 
N 



3. The number of terms in Eq. (|IBD and Eq. (|19j) are both 2 



M 



After factoring out the integrals over the momenta and positions of the M — N solution layer 
particles, Eq. (O) becomes 



Z tot = E^ =0 ^ AM (M!)- 1 (/d d rj M ^(/dVexp(-/3p 2 /2m)) y 
x / d d pi • • • d d p N d d n ■ ■ • d d rWexp (-(3 [H N + $(JV, M, fi, (3, V)}) . 
This equation can be easily rewritten, 



M-N ( M 



(20) 



2 



tot 



^M 



X 

£ 



{V) M-N ( 



2TTmk B T \ d ( M - N ">/ 2 e X p(-f3<S>(N,M,ti,l3,V)) 
h 2 ) (M-N)l 



h -Nd( N iyi j d dfo . . . d^dVx • ■ • d d f N exp (-PH N 

(M-N)] 



(21) 



M 

N=0 



<yf- N 

Equating Eq. fl2T|) and Eq. (|T5|) now yields 

V \ M ~ N exp (-P$(N, M, /x, /3, 7)) 



N 



A d 



where 



A 



(M — N) ! 



/i 2 



exp (/3/ziV) 



1/2 



(22) 



(23) 



is the thermal wavelength. 

The value of the correction factor is thus 

$(A, M, /I, p, V) = -fiN + k B T 



(M-N)\n(J^-ln((M-N)\) 



(24) 



This correction term may be included directly into the Hamiltonian, or we may redefine the 
zeros of the energy and chemical potential, 



Hi 



tot 



Htnt - k R TM\n 



(a«0 



(25) 
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and 

fi' = ti + k B Thi(^j . (26) 
This yields the effective Hamiltonian 

M 2 M-l M M 

Kot = ESr + E E U int(\n - r 3 \)c t c 3 + J2 UU?i)c* - l/N - k B T\n((M - N)\) . (27) 
i=i zm i=i j=j+i i=i 

Monte Carlo energy differences are then calculated with this effective Hamiltonian. There is 
no need to treat the particle momenta in a system with no momentum-dependent forces, and 
we do not change the kinetic energy during the Monte Carlo procedure. 
The Monte Carlo proceeds in the following way: 

1. Choose one of the M particles at random. Call the chosen particle i. 

2. Randomly choose to diffuse or adsorb/desorb particle i. The relative probability of choos- 
ing diffusion or adsorption/desorption is not particularly important, and for simplicity, 
the probability of choosing diffusion is set to 1/2. 

• If diffusion is chosen, choose the new r. \ 1 = f, \ + dr, where ' is the new position 
and dr is a uniformly chosen displacement within a d- dimensional sphere of radius 
R. 

• If adsorption/desorption is chosen, change the layer occupation variable, q, such 
that if Cj = 0, then d i = 1, and if c, = 1, then d { = 0. 

3. Calculate the energy change using Eq. (p7|) . 

4. Accept the new configuration with any detailed-balance satisfying rate, like the Metropolis 
rate. 

5. Go to 1. 

The energy differences are given here for the four different cases. 

(i) q = c'j = 1, the particle diffuses within the surface layer. 

M 

A?4t = E [ c i {Uint(\n ' - r 3 \) - U- mt {\n - r 3 \))} + U cxt (r t ') - U cxt (r t ) (28) 

(ii) Cj = c'j = 0, the particle diffuses within the solution layer. 

A?4t = (29) 
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(iii) Q = 1, c- = 0, the particle desorbs. 

M 

AH' tot = -C/ ex t(n) + // - k B Tln(M -N + l)-^ c 3 U hlt (\n - (30) 

(iv) Q = 0, = 1, the particle adsorbs. 

M 

AH[ ot = C/ e xt(n) - // + A; B Tln(M - N) + ]T c^ int (|rl - r ^ | ) (31) 

Here, TV represents the initial number of surface-layer particles, corresponding to the configu- 
ration at step (1). 
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Table 1: DFT calculations of the corrugation potential for Br/Ag(100) are shown in column 4. 
More negative values indicate stronger binding, and (x, y) = (0, 0) indicates the on-top site. 
Points related by symmetry are not listed. Column 5 lists the difference between the DFT 
calculations and the sinusoidal approximation (A = 150 meV) of Eq. (|I|). 



x/a 


y/a 


site 


U D FT(x,y)/meV 


[U DFT (x,y) -U(x,y)]/meV 








on-top 


305 


155 





0.25 




206 


131 





0.5 


bridge 








0.25 


0.25 




121 


121 


0.25 


0.5 




-54 


21 


0.5 


0.5 


hollow 


-154 


-4 



Table 2: Effective Arrhenius barriers for different corrugation amplitudes and damping pa- 
rameters. The errors are estimated from the linear fits shown in Fig. ^. 



A [meV] 


Effective Barrier [meV] 


a = 1 


a = 10 


a = 100 


100 


85.2 ±0.8 


100 ±1 


98 ±2 


150 


125.0 ±0.7 


148 ± 1 


144.9 ± 0.8 


200 


164 ± 1 


198 ± 1 


199 ±1 



Table 3: Effective attempt frequency for different corrugation amplitudes and damping pa- 
rameters. The errors are estimated from the linear fits shown in Fig. ^| The time scale for one 
attempt is u^ 1 . 



A [meV] 


Effective v 


s- 1 ] 


a = 1 


a = 10 


a = 100 


100 


1.14 ±0.04 x 10 12 


2.3 ±0.1 x 10 12 


4.4 ±0.5 x 10 11 


150 


1.31 ±0.04 x 10 12 


2.9 ±0.2 x 10 12 


5.8 ±0.2 x 10 11 


200 


1.52 ±0.08 x 10 12 


3.5 ±0.1 x 10 12 


9.1 ±0.6 x 10 11 



19 



Fig. 1 The corrugation potential for Br adsorbed on Ag(100), approximated as a two- 
dimensional sinusoidal function. The grayscale indicates the Br binding energy, 
with darker shades indicating more favorable binding and lighter shades indicating 
less favorable binding. The Ag surface atoms giving rise to the corrugation potential 
are shown as circles with diameter a = 2.889 A [0. Some of the circles have been 
removed for clarity. Top, bridge, and hollow sites are labeled. 

Fig. 2 Visualization showing both the solution layer and the surface layer. The spheres 
indicate Br ions, where the smaller inner sphere is shown for graphical purposes 
to indicate the center of the ion, and the outer wireframe shell has radius i?; on = 
1.94 A |fL6| , corresponding to the ionic radius of Br~. Br in the solution are confined 
to the upper plane and are non-interacting. Br on the surface are confined to the 
lower plane and interact with both the substrate corrugation potential and with 
other adsorbed Br particles. The magnitude of the corrugation potential is shown 
by the grayscale shading on the lower plane. Darker shades indicate more favorable 
locations. 

Fig. 3 Lateral interaction energy as a function of the lateral separation between a Br pair. 

See Eq. (||). The vertical lines indicate important distances as indicated in the figure. 

Fig. 4 Root mean square (rms) displacement of Br adparticles from the four-fold-hollow 
sites vs coverage 0, based on an off-lattice MC simulation in the canonical ensemble. 
The system size is L = 32, and the results for three different values of A are shown. 
The displacement is largest for low coverages and small A and smallest for high 
coverages and large A. 

Fig. 5 Off-lattice grand-canonical MC isotherms for L = 32 with three different values of 
A. The differences are most pronounced in the disordered region around ~ 1/4. 

Fig. 6 Single-particle trajectory from a Langevin simulation with a = 10 for A = 150 meV 
at room temperature. The gray-scale indicates the corrugation potential, the heavy 
curve indicates a simulated particle trajectory, and the dashed lines indicate the 
unit-cell boundaries used in defining crossings. 
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Fig. 7 The crossing-time probability density from a Langevin simulation with a = 10 and 
A = 150 meV for single-particle crossings at room temperature. The linear region 
in the plot indicates the large-r Arrhenius behavior. Significant deviations from this 
behavior are seen at small r. 

Fig. 8 Directional crossing probabilities from a Langevin simulation with A = 150 meV 
and three different values of a for single-particle crossings at room temperature. 
For small r, the backward crossings indicate crossing/recrossing events, and the 
forward crossings at small r indicate running events. The perpendicular direction 
includes both left and right crossings and has been divided by two for comparison 
with the other crossing probabilities. The figures show (a) a = 1, (b) a = 10, and 
(c) a = 100. 

Fig. 9 Plot to determine the Arrhenius parameters from the Langevin single-particle cross- 
ing times. The slope gives -Ebam and the intercept gives u, which are listed in Table |] 
and Table [3], respectively. The figures show (a) a = 1, (b) a = 10, and (c) a = 100. 
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